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The cascade-shell model of turbulence with six real variables originated by Gledzer is 
studied numerically using Mathematica 5.1. Periodic, doubly-periodic and chaotic solutions 
and the routes to chaos via both frequency-locking and period-doubling are found by the 
Poincare plot of the first mode v%. The circle map on the torus is well approximated by 

(N ■ 

the summation of several sinusoidal functions. The dependence of the rotation number on 
the viscosity parameter is in accordance with that of the sine-circle map. The complicated 
i bifurcation structure and the revival of a stable periodic solution at the smaller viscosity 

parameter in the present model indicates that the turbulent state may be very sensitive to 
the Reynolds number. 
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£^ . Many of turbulent phenomena in motions of fluids are considered to be understandable 
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(NS) equations. Recent progress of computer hardware and software, however, may still not 
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be enough to make clear the route and the characteristics of turbulence. Another prompt way 
to smatter turbulence is to omit the theoretical derivation of nonlinear models, to give up the 
full numerical simulation and to fall back on simpler models analogous to the NS equation. 
In this context, the cascade-shell models by Gledzer (1973) 1 and its complex version (Ohk- 



itani and Yamada (OY 1989), 2 so-called the GOY model; see also Frisch (1995), 3 Kato and 
Yamada (2003), 4 Biferale(2003) 5 for further references) have been studied numerically from 
the viewpoint of turbulence statistics. In the present investigation, the bifurcation approach 
is made to the Gledzer's cascade-shell model of turbulence with six real variables. Biferale et 
al. (1995) 6 has found the limit-cycle and torus attractors in the GOY model with a fixed vis- 
cosity along with loss of the stability of the Kolmogorov 1941 fixed point when the parameter 
e related to the helicity exceeds critical values about 0.386 and 0.396, respectively. In contrast 
with the bifurcation found by Biferale et al. (1995) 6 in GOY interpreted as the Ruelle-Takens 
scenario, a detailed study of the circle map in the Gledzer model in this paper indicates that 
the bifurcation of torus attractors is in accordance with that of the sine-circle map. Another 
similar study of the five-mode truncation model of the NS equation is made by Franceschini 
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and Tebaldi (FT 1979), 7 where the period-doubling and the symmetry-breaking bifurcations 
are found but the quasi-periodic route is out of their result. 

The Gledzer's original model of real variables Vi(t) for i = 1, • • ■ , n is given by: 

dvi 

— = Vi = CliV i+ lV i+ 2 + C 2 ,iVi-lV i+ l 

at 

+c 3ii fi_ii;i__2 - vkfvi + fi, (1) 

where Vi (fi) is the velocity (forcing) of the i-th mode in the space of a discretized wavenumber 
fej. v denotes the kinematic viscosity, the forcing is assumed to be time independent, and 
vq = V-\ = v n+ i = v n+ 2 = 0. In order to conserve the energy E = Y^i=x v 1 an d the enstrophy 
Y17=l ki v i m t ne case °f v = fi = as Gledzer (1973) required for the model to be analogous 
to the two-dimensional turbulence, the coefficients of the nonlinear terms c,-j for j = 1,2,3 
and i = 1, ■ • • , n need to satisfy the following relations: 

1,2 _ L2 

C 3 ,i+2 - -T2 _ .2 C M- V^i 

^i+l %+2 

If the wavenumber and the nonlinear coefficient en are selected as ki = en = koq l , the above 
relations become 

c 2 ,i = -/0fe_i, (4) 



C3,i = (/9 - l)fci_ 2) (5) 

where 

/3 = l + g- 2 . (6) 

However, we switch the parameters chosen by OY, (5 = 1/2 and q = 2 in iflj). fcj and ci^, 
which do not satisfy @, model ling the 3D turbulence in the sense that the conservation law 
holds only for energy, not for enstrophy. 

Assuming that the forces are applied only on the first and second modes fx = fi = 1 and 
fi = for i > 3, the equation for the total energy with the forcing and the viscosity becomes 

n 

E/2 = f lV x + f 2 V2-^kfvf. (7) 

8=1 

If we consider the sufficiently large values of such that \vi\ > \fi\/{vkf) for the forced modes 
i = 1,2, the right hand side becomes negative and the solution is proved to be finite. The 
present model deals with real variables and therefore it is possible to reduce the dimensions 
of the model to the half with the same extent of the wavenumber, compared with the GOY 
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model. 



Letting ko 



2 4 , the model becomes 



Vi + \V i+ 2 ~ 




(8) 



The present model possesses a single parameter u, changed between the range about 5-10 2 - 
5-10" 3 . Three ways of giving the initial condition are considered in order to examine multiple 
stable states. Case (I) is the origin; i>i(0) = 0. Case (II) [or (III)] is the final computed values 
of the slightly smaller [or larger] value of the viscosity v. 

The Linux version of Mathematica 5.1 is used for the main part of the numerical integra- 
tion. According to the manual of Mathematica, the general-purpose ODE solver developed by 
Hindmarsh (1983) 8 is adopted in the NDSolve command. The numerical solution is constructed 
checking its local convergence up to the machine precision of order 10 -16 . The weakest point 
is that it consumes large memory. The number of the modes is fixed as n = 6 in the main part 
of this study. The computer has 1GB memory, its CPU is AMD Athlon XP 2000+ and its 
operating system is Vine Linux 3.2. A Fortran program with the fourth-order Runge-Kutta 
scheme with a fixed time step is also coded to check the bifurcation diagram in Figure 2 and 
no inconsistency is found between results by Mathematica and Fortran. 

The number N(n) of fixed points (FPs) of (JHJ) including complex numbers are computed 
by Mathematica and shown as JV(4) = 5, JV(5) = 9, JV(6) = 25, JV(7) = 29, N(8) = 61, 
N(9) = 129 and iV(10) = 177, but the number of real fixed points is only 5 for n = 4, 5 and 
3 for 6 < n < 10 at v = 10~ 2 . 

Figure 1 shows the bifurcation structure of the FPs; the stable (unstable) FPs are denoted 
by thick (thin) curves. The Mathematica program is coded such that the number of modes 
is arbitrary. The fact that the bifurcation structure is quantitatively the same between n = 6 
and 8 suggests the sufficiency of the six-dimensional system in the region 0.005 < v < 0.05. 
The number of stable FPs is 1 for v > 0.03983(3), 0.03482 > v > 0.02059(5), 0.01287 > v > 
0.01202(5), 2 for 0.02058 > v > 0.01288(5) and for 0.03982 > v > 0.03483(3), 0.01201 > v. 
The number of total FPs is shown in the above parenthesis; 3 for 0.05 > v > 0.03483, 
0.01167 > v > 0.00967, 5 for 0.03482 > v > 0.01168, 0.00968 > v > 0.00791, 7 for 0.00790 > 
v > 0.00355 and 9 for 0.00354 > v. There is a tendency that as v decreases, the number of 
real FPs increases. 

Figure 2(a) shows the bifurcation diagram of attractors by plotting values of v\ at the 
local maxima of v\ (v\ = and v\ < 0) between t = 9 ■ 10 3 and 10 4 , actually computed by the 
interpolation of three adjacent points with the time interval dt = 0.01. The initial condition 
is Case (I). The curve in the center of Fig. 2(a) is therefore due to the stable FP. 

The parameter region including the stable doubly-periodic state is examined in detail using 
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the initial condition of Case (II) and (III). The stable doubly-periodic solutions are confirmed 
to be generated through the supercritical Hopf bifurcation at about v = 0.03834 as shown 
in the enlarged Figure 2(b) and Figure 2(c). The doubly-periodic solutions coexist with the 
periodic solutions with period 6T, on which the incomplete period-doubling bifurcation occurs 
as seen in Fig. 2(b). Here, the period mT means that the periodic solutions give m points on 
the present Poincare plot. The torus breakdowns at about v = 0.03785 where the 6T-periodic 
solution is still stable. Fig. 2(c) shows that there are windows of frequency-locking, which 
is one of the remarkable features of the sine-circle map. According to the current resolution 
in Figure 2(d), the 6T-periodic solution changes chaotic suddenly at about v = 0.03724, 
indicating the intermittency route. 

The left side of the chaotic parameter region in Fig. 2(a) is enlarged in Figure 2(e). 
Periodic windows are also observed among the chaotic solutions. A sequence of period-doubling 
bifurcations (Feigenbaum 1978) 9 of orbits is also observed in Figure 2(f), although it is in fact 
a revival of the periodic solution from chaos since the horizontal axis denotes the viscosity 
and the inviscid limit v — > turns to the left in the Fig. 2(f). Numerically computed examples 
of attractors projected on the (vi,V2,vz) space for various values of v are shown in Figure 3. 

In order to examine the circle map of the quasi-periodic solution, the pair (v^ , {^™ +1 ^) is 
considered, where = — < v[ n ^ >, is the n-th plotted value of v\ and < • > denotes 
the time average. The points lie on a closed curve in the Poincare section, showing the evidence 
of the doubly-periodic torus motion in the original 6D space. Denoting the 27r-normalized 
argument of v[ n) + ivf +1) by 9 n (0 < 9 n < 1), the one-dimensional map 9 n +i = f(9 n ) can be 
constructed. 

In many studies of the circle map, the sinusoidal function is considered, but strictly speak- 
ing, the function is likely to deviate from it in real chaos systems of ODEs. We seek an 
approximation of the circle map of the form 



by applying the FindFit command of Mathematica to the numerically obtained data of the 
Poincare plot with M = 10. Figures 4(a-c) show the dependence of Kj(j = 1, • • • , M) and 
£1 on the viscosity parameter. We observe the nonvanishing values of K2, K4, Kq even at the 
onset of the stable torus; i.e. the second Hopf bifurcation. 

The numerical computation of the rotation number p on the torus is also possible by 
counting the number JV^ of the decrease of 9 n due to the modulus 1 in Eq. Q. A naive 
numerical approximation is p = N^/N, where N is the number of the total Poicare plot, 
denoted by small dots in Figure 5. A more elaborate method uses the power spectrum of 



9~n+l — f(@n) — 9 n + Q 




M 



mod 1 



(9) 
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vi(t). Choosing the suitable pair of two frequencies giving the maximum of the spectrum, we 
can identify the rotation number, denoted by larger dots in Fig. 5. These two results are in 
agreement within the resolution of the numerical computation. The qualitative similarity is 
observed with that of the sine-circle map, showing the monotonic decrease of p and the flat 
region, i.e. the evidence of the incomplete devil's staircase. 10 Since our system has only one 
parameter, we just need to imagine the situation that two parameters (0, K) in the sine-circle 
map are changing simultaneously, as shown in Figures 4. 

In summary, the two different routes to chaos coexist in the single Gledzer's model of 
six real variables. As noted by E. Ott 11 on page 199, the quasi-periodic route to chaos does 
not necessarily imply the bifurcation from double periodicity directly to triple periodicity or 
chaos, so-called the Ruelle-Takens scenario. The study of the sine-circle map indicates that 
frequency-locking leads to periodic solutions, which may become chaos through the period- 
doubling route or the other ones like symmetry-breaking or intermittency. The author has 
tested the Langford 12 equation, well known to have torus attractors, and finds a similar 
scenario suggested by Ott. 11 

According to Figures 4 and 12 in OY, 2 10 positive Lyapunov exponents have been found 
for the 48-dimensional model of 24 complex variables and the dimension of the corresponding 
attractor should be very large. Their parameters v = 10 -8 and / = (1 + i) x 5 x 10 -3 
correspond to the normalized viscosity v' = u/Kef = 2 x 10~ 6 in our model with / = 1. Their 
computation time t = 400 leads to the normalized time t' = (Kef)t = 2 in our case with n = 6. 
The computation of Lyapunov exponents of our model, as well as motions of eigenvalues of the 
Jacobian, have been made but is not shown in this Letter. The triply-periodic motion has not 
yet identified in our model. It is still open whether the low-dimensional chaos shown in this 
study is connected or not with the state of turbulence which is considered to be independent 
of the parameter z/, as we increase the number of n and decrease v. From the viewpoint of the 
bifurcation, it requires careful analysis but the gap will be filled as the power of the computer 
increases. 

The author is grateful to Prof. T. Yamagata for encouragement, to Prof. H. Sakai for 
permission to use the file server, to Prof. W. F. Langford for sending me the related reprints 
on the Langford equation, and to Prof. L. Biferale for letting the author know his work before 
proofreading. 

Figure 1. Bifurcation diagram of the fixed points (FPs) between 5 • 10 -2 > v > 5 • 10~ 3 . 
Stable (unstable) FPs are shown by thick (thin) curves. 

Figure 2. Bifurcation diagram of attractors shown by values v i at the local maxima of v\ 
{v\ = and v\ < 0). 
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(a) The range of the parameter v around 5 • 10~ 2 — 5 • 10~ 3 and the initial condition is the 
origin. 

(b) The parameter is 0.0372 < v < 0.0409 and the initial condition is Case(II) and (III) 
such that the multiple stale states of the double periodicity and the 6T-periodicity are clearly 
captured. 20000 points are selected randomly for minute drawing. 

(c) The enlarged Figure of (b) at 0.0378 < v < 0.03834 showing the second supercritical 
Hopf bifurcation, the doubly-periodic and frequency-locked periodic solutions (especially, 11:9, 
17:14, 23:19 and 29:24 resonance), in addition to another stable 6T-periodic solution. 

(d) The parameter is 0.0372 < v < 0.0384. The chaotic solution emerges suddenly at 
about v = 0.03724, indicating the intermittency route. 

(e) Chaotic attractors with periodic windows at 0.006 < v < 0.012. 

(f) The enlarged figure of (e) at 0.006734 < v < 0.00674, showing the evidence of the 
period-doubling bifurcation. 

Figure 3. Projections of attractors on the (vi, V2, V3) space at v = (a) 0.0395, (b) 0.0379, 
(c) 0.0378, (d) 0.037, (e) 0.0355, (f) 0.0118, (g) 0.0103, (h) 0.00883, (i) 0.0076, (j) 0.0074, 
(k) 0.0064 and (1) 0.006. (a) and (1) are the periodic solution of the period T, (k) 2T, (f) 
3T, (h) 4T and (g), (c) 6T. (b) is the doubly periodic and (d), (j) are the chaotic solutions, 
respectively, (e) appears to be close to the heteroclinic orbit. 

Figure 4. The dependence of the parameters in the circle map K{ with even i (a), odd i 
(b) and O (c) on the viscosity v. 

Figure 5. The numerically obtained dependence of the rotation number p on the viscosity 
v. The small (large) dots denote p by the method counting the number of decrease of 9 n due 
to the modulus (computing the ratio of two frequencies giving the maximum of the power 
spectra of v 1). 



6/EH1 



J. Phys. Soc. Jpn. Letter 

References 

1) E. B. Gledzer: Sov. Phys. -Dokl. 18 (1973) 216. 

2) K. Ohkitani and M. Yamada: Prog. Theor. Phys. 81 (1989) 329. 

3) U. Frisch: Turbulence (Cambridge University Press 1995) 

4) S. Kato and M. Yamada: Phys. Rev. E 68 (2003) 025302. 

5) L. Biferale: Ann. Rev. Fluid Mech. 35 (2003) 441. 

6) L. Biferale, A. Lambert, R. Lima, and G. Paladin: Physica D 80 (1995) 105. 

7) V. Franceschini and C. Tebaldi: J. Stat. Phys. 21 (1979) 707. 

8) A. C. Hindmarsh: Scientific Computing 8 (1983) 55. 

9) M. J. Feigenbaum: J. Stat. Phys. 19 (1978) 25. 

10) M. H. Jensen, P. Bak and T. Bohr: Phys. Rev. A 30 (1984) 1960. 

11) E. Ott: Chaos in Dynamical Systems (Cambridge University Press 1993) 199. 

12) W. F. Langford: Nonlinear Dynamics and Turbulence (Pitman Advanced Publishing Program 
1983) 215. 



7/DJ 



200 



150 



100 



50 








13 



en 
o 
o 



P 



0.01 



0.02 



0.03 



0.04 



0.05 



t- 1 



13 



en 
o 
o 



13 




10 



15 Vi 

era 
(3 

n 

o 

to 








F2V 



: ilL 



0.005 



0.01 



0.015 



0.02 



0.025 



0.03 



0.035 



0.04 



t- 1 



Phys. Soc. Jpn. 



Letter 




Figure 2b 



lO/HSl 



Phys. Soc. Jpn. 



Letter 




Figure 2c 




Figure 2d 



12/Trgi 



J. Phys. Soc. Jpn. 



Letter 




I 



I 



Figure 2e 



13frai 



Phys. Soc. Jpn. 



Letter 



MX 



I r 



_i i_i 



>,l, -, 



i n, 1 ■ ■ ■ 

llHii fl 

if i iii ii 

ii j iiii ij 

ii i tit} it 

iii iiii ii 

iii iiii ii 

iii ill ii 
iii iii ii 

i MM ■ 

■ i ii i 

ii iii 

i \ ii i 

■ MM 

j i ii i 

j j jjj 

■ i iii 

■ i iii 
i i iii 

i i iii 

j j jfj 
i y i j 

iiii 
iiii 
iiii 
iiii 

iiii 

i i ii 

i i ii 



III 

iii 



1 



•n 

s 

d 



§3 



00 

m 

& 
O 
O 

d 



O 



5 



Figure 2f 



14/Q33 



J. Phys. Soc. Jpn. 



Letter 




Figure 3 
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